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Abstract 

Polyconvexity is one of the known conditions which guarantee existence of so- 
lutions of boundary value problems in finite elasticity. In this work we propose a 
framework for development of polyconvex strain energy functions for hyperelastic 
materials with cubic anisotropy. The anisotropy is captured by a single fourth order 
structural tensor for which the minimal polynomial basis is identified and used for 
the formulation of the strain energy functions. After proving the polyconvexity of 
some polynomial terms, we proceed to propose a model based on a simple strain 
energy function. We investigate the behavior of the model analytically in one di- 
mension and numerically in two and three dimensions. These investigations allow 
us conclude that the model possesses the physically relevant directional properties, 
in particular, strains in different directions are different and the lack of any loading 
symmetries leads to development of "shear" stresses and displacements. 
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1 Introduction 

Anisotropic engineering materials exhibit directionality in their mechanical 
characteristics even when subjected to very large strains. These materials rep- 
resent a wide range of applications in composites and crystals as well as in 
bio-mechanical systems. Although some advances have been made in the direc- 
tion of characterizing simple cases of anisotropic material behavior through the 
proposal of phenomenological models respecting the applicable mathematical 
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theories, the field is far from completion. In general most of the phenomenolog- 
ical studies lack detailed analyses on the mathematical properties of the pro- 
posed models. Even in [1], where an analysis of general convexity conditions 
for transversely isotropic materials is extensively presented, the relationship 
between the numerous proposed functions with their physical counterparts is 
still to be fully developed. 

The material symmetries of an oriented continuum impose definite restrictions 
on the form of constitutive relations. The procedure used for the construction 
of constitutive models must from the very beginning assure that the equations 
are written in a proper manner which refiects the material symmetries. Fur- 
thermore, the final goal of the procedure is the development of a mathematical 
framework that satisfies conditions guaranteeing the existence of solutions for 
models that lack the standard regularity properties assuring existence and 
uniqueness. Indeed, uniqueness should not be required because it precludes 
description of some important physical effects of great interest, for example, 
buckling (in this respect see [2]). The procedure, presented in this work, follows 
the approach laid out in [1]. 

The fundamental aim from a mathematical perspective is to guarantee the 
existence of solutions. Existence of minimizers of some variational principles 
in finite elasticity is based on the concept of quasiconvexity (introduced by 
Morrey in [3]), which ensures that the functional to be minimized is weakly 
lower semi-continuous. Unfortunately, quasiconvexity gives rise to an integral 
inequality, which is extremely difficult to handle due to its global character. 
Therefore we turn to the more practical concept of polyconvexity ( [2] ) which 
can be verified locally. 

The increased complexity of observed mechanical behavior of anisotropic ma- 
terials requires invariant formulations of anisotropic constitutive laws. The 
theory of tensor function representations constitutes a rational procedure for 
consistent mathematical modeling of complex anisotropic material behavior. 
A particularly strong push in that direction is the work [4] by Weiss who 
introduced an exponential function in terms of the mixed invariants. 

As extensively presented in [1] , the complex mechanical behavior of elastic ma- 
terials with an oriented internal structure at large strains can be described with 
tensor-valued functions of several tensor variables: the deformation gradient 
and a few additional structural tensors. The follow-up strategy is to construct 
constitutive models with an invariant form of the strain energy function. The 
general forms of tensor-valued functions have been derived in the form of of 
representation theorems for tensor functions ([5]). The type and minimal num- 
ber of the scalar variables entering the constitutive equations are also known. 
The interested reader should consult [6,7,8,9,10,11] for details. 
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In this paper we are concerned with the problem of determining the general 
form of scalar-valued polynomial expressions, and for that reason we make 
use of the concept of integrity basis. So far in the literature only the simplest 
form of material anisotropy represented by transversely isotropic materials 
with a single preferred principal direction, has been extensively developed 
following these lines. In this work we develop a procedure for the construction 
of polyconvex free energy functions for cubic crystal systems. These cubic 
crystal systems present three orthogonal principal directions, giving rise to a 
considerable increase in the complication of the mathematical machinery to 
be dealt with. The main difference with previous works comes from the need 
to use a fourth order structural tensor to characterize the material symmetry 
group. The need for the fourth order tensor comes from our desire to use a 
single structural tensor. We will make use of results obtained by Zheng ([12]) 
on the single structural tensors characterizing the crystal classes. 

To summarize, this work presents a large deformation mathematical model for 
anisotropic materials with cubic symmetry. 

The paper is organized as follows. In Section 2 we present the basic notation 
and and review some kinematics relations at finite strains to be used in the se- 
quel. After that we focus on the presentation of the mathematical framework 
for hyperelastic materials which guarantee a priori some meaningful physi- 
cal conditions, in particular, the material frame indifference and the material 
symmetry conditions; it will be shown that these two conditions require the 
introduction of the concept of structural tensors. Section 3 is concerned with 
the application of the concepts of hyperelasticity and structural tensors to 
the particular case of a material formed by cubic crystal. After characteriz- 
ing the material symmetries associated with the cubic crystal anisotropy by 
means of the appropriate structural tensor, we present a procedure to build 
up free energy functions for cubic crystal materials that fulfill the appropri- 
ate mathematical requirements, more specifically the polyconvexity condition. 
The proposed functions have the invariants of the deformation gradient and 
the structural tensor as arguments. This approach requires the concept of 
polynomial integrity basis, which is also presented. The representation for the 
stresses and the tangent matrix is given in detail. A model fulfilling all the 
requirements mentioned above is finally proposed in section 4. Two conditions 
are added to fully determine the problem and relate it to the physical data. 
These conditions are the stress-free reference configuration and the linearized 
behavior near the natural state C = 1; these are dealt with in section 5. We 
consider the behavior of the proposed model in one dimension in the next 
section, where its physically desirable stress-strain response can be fully ap- 
preciated. A short summary of the variational and finite element formulation 
is given in section 7. Section 8 presents numerical results obtained from sim- 
ulation examples using the proposed model. Finally, in section 9 the main 
conclusions of the present work are summarized. A few appendices have been 
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added at the end to encompass some of the derivations. 



2 Foundations of continuum mechanics 

In the following we consider the class of hyperelastic materials for which we 
postulate the existence of a free energy function. The resulting constitutive 
equations must fulfill some requirements that naturally arise from physical 
considerations of response invariance of the material under arbitrary coordi- 
nate system transformations. It will be shown that the requirement on the 
constitutive functions of anisotropic solids to satisfy the material frame indif- 
ference will force these functions to be isotropic tensor functions. Therefore 
the material symmetry condition cannot be accomplished simultaneously with 
the material frame indifference. In order to resolve this issue we will make use 
of the concept of structural tensors. Structural tensors increase the number of 
arguments of the energy functions, enabling the model to account separately 
for the material symmetries. 

2.1 Notation and kinematics 

In this section we briefly present the notation and main results corresponding 
to kinematics in the standard theory of continuum mechanics. The theory 
presented here is based on a material formulation. 

The movement of a continuum body can be seen as a family of configurations 
ordered by the time parameter. Thus, for every t e [0, T] C SfJ"*", the application 
0( : i? — s> 5 C 3?^ is a deformation which transforms the reference configuration 
B in the configuration S at time t. Then, x = 0t(X) = 0(X, t) identifies the 
position X of point X at time t. We will follow a lagrangian description of 
the motion, which implies that the material coordinates of a point, {^a}? are 
taken as independent variables. It is usually called material description of the 
motion. The deformation gradient F is defined as 

F = V0t(X) 

with the jacobian J = det(F) > 0, defined positive as a condition to prevent 
material interpenetration. 

Let F denotes the material time derivative of the deformation gradient. It is 
identical to the material velocity gradient, i.e. 
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The deformation gradient F can be used to form the right Cauchy-Green 
tensor, which corresponds to the chosen strain measure, i.e. 

C = F^F. (1) 

In general spaces all inner products appearing in the former and latter deriva- 
tions should be properly constructed taking into account the corresponding 
space metric, defined as a symmetric second order covariant tensor, and de- 
noted by G for the reference configuration, and by g for the deformed config- 
uration. In this paper we will restrict ourselves to the case of euclidean space 
with cartesian coordinates in which the metric tensors become the identity, 
and therefore they will not explicitly appear in the calculations. 

2.2 Hyperelasticity and invariance conditions 

As mentioned in the introduction, we will focus our study on hyperelastic 
materials. They are an elastic materials class which postulates the existence 
of a stored free energy function ^/;(X,F, ■). The energy function ip depends 
on the point in the reference configuration X, the deformation gradient F, 
and an additional tensor, which characterizes the material anisotropy. We will 
restrict ourselves to perfectly elastic materials, which means that the inter- 
nal dissipation Dint is zero for every admissible process (see [13]). Following 
a standard argument, the constitutive equations relating the stresses to the 
energy function are obtained by evaluation of the Clausius-Planck inequality 

Ant = P : F - = (P - 9f^) : F > ^ P = 9f7/; (2) 

where the thermal effects have been neglected and P is the first Piola-Kirchhoff 
stress tensor. 

The principle of material frame indifference requires the invariance of the 
constitutive equation under rigid body motions superimposed onto the current 
configuration, i.e., under the mapping x Qx the condition ?/'(F) = ?/;(QF) 
holds for every Q in the special orthogonal group SO (3). 

As shown in [14] the requirement that the constitutive equations fulfill the 
principle of material objectivity yields the functional dependence ip = ip{C) = 
ip{C{F)), i.e., every dependency on the gradient F can be properly substituted 
for a dependency on C. Now, considering the relation between the first and 
second Piola-Kirchhoff stress tensors, together with the dependency of the 
energy function on C and expression (2), we can obtain the relation between 
stresses measured by the second Piola-Kirchhoff stress tensor and the energy 
function. Thus, we have 

S = F^ip = F~1(9f^ = F-1((9c^(9fC) (3) 
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and considering the symmetry of C, we deduce that dc'ipdFC = 2Fdcip, which 
carried into the expression for S above gives 



S = 2dc^. 



(4) 



The anisotropy of a material can be characterized by the material symmetry 
group Gm with defined respect to a local reference configuration. The elements 
of Gm are those transformations Q that give an invariant material response, 
i.e., they are superimposed rotations and reflections on the reference config- 
uration which do not influence the behavior of the anisotropic material, thus 



The conditions (5) establish that the function tp and the tensor P are Gjv/— in- 
variant. In general we have Gm C S0{3), so the material symmetry group 
corresponds to a subgroup of the whole special orthogonal group 5*0(3), and 
only in the case of an isotropic material both groups coincide. This last fact 
gives rise to a problem regarding two conflicting requirements: from one side 
the functions should be invariant only under transformation belonging to 
Gm, reflecting the material anisotropy, and from another side their formu- 
lation should be transformation independent, so that the representation is 
coordinate- free, i.e., the Q's in (5) should belong to 5*0(3). To summarize, 
the material symmetry condition requires the use of an isotropic function, but 
at the same time that leads to loss of the information concerning the material 
anisotropy. 

It has to be emphasized, however, that so far we have been considering only 
constitutive functions dependent on one argument, the tensor C, and this 
points out to one possible approach to meet both requirements. It will be 
shown that both requirements can be satisfied simultaneously by extending 
the tensorial argument list of the energy functions, thus obtaining an isotropic 
function embodying the anisotropy information. This approach is put into 
practise by means of structural tensors. 

2.3 Isotropic tensor functions for anisotropic material response. Structural 



As it has been shown in the previous section, the constitutive equations can 
be deduced from a free energy function, but we faced the problem of charac- 
terizing anisotropic materials with a dependency on C only. We saw that it 
was not possible to have both the anisotropy and the invariance under any 
spatial rotation and reflexion captured in that manner. The idea behind the 
structural tensors is to be able to have an isotropic tensor function, i.e., one 



^(CQ) = V^(C) VQgGm 
P(CQ) = P(C)Q VQ e Gm. 



(5) 



tensors 
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which is invariant under any rotation in space, but at the same time to retain 
the symmetry information characterizing the anisotropy of the materiaL Both 
conditions of rotation invariance and anisotropy can be properly fulfilled by 
adding more tensors as additional arguments in the free energy function. 

The structural tensors are useful in obtaining irreducible and coordinate-free 
representations for anisotropic tensor functions because they characterize the 
spatial symmetry group. The concept of structural tensors was introduced by 
Boehler in [7]. The characterization of the symmetry group is done in the 
following sense: the tensors ^, C are said to be the structural tensors of the 
spatial symmetry group G m if and only if 



Basically, the effect of the structural tensors can be captured in the difference 
between the following two statements 

• The relation ijj{C,C,) = ^/'(Q'^CQ, Q^.^) for the free energy function, and 
consequently for the corresponding stress tensor, Q'^S(C, ,^)Q = S(Q'^CQ, Q*^), 
holds for VQ G SO (3), which means that the function is an isotropic scalar- 
valued tensor function. 

• The relation ip{C,^) = 'i/'(Q'^CQ,^) for the free energy function, and the 
corresponding stress tensor, Q'^S(C,^)Q = S(Q'^CQ,^), holds for VQ G 
Gm, which means that the function is an anisotropic scalar- valued tensor 
function. 

On a more intuitive level, the difference between these two statements is as 
follows: in the first one the deformation and the "body structure" are both 
rotated with Q, while in the second statement only the deformation is rotated 
because the structural tensor C, has the appropriate symmetry to rotations 
with Q G Gm- 



3 Free energy function for cubic materials 

As we showed in the previous section the final aim in the proposal of a model 
is to construct energy functions invariant under the appropriate symmetry 
groups reflecting the underlying material anisotropy. A direct way to do that 
is by means of functions dependent on the invariants of the right Cauchy- 
Green tensor and structural tensors, which ensures that the functions to be 
constructed are also invariant under the proper symmetry group, retaining in 
this manner the material symmetries of the body of interest. In particular we 
are interested in proposing polynomial type energy functions. In order to keep 
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the complexity to minimum we will make use of the minimal set of independent 
invariants of the deformation and structural tensors. This minimal set is called 
a polynomial basis and further details about it can be found in [8]. 

In the subsections to come we first present the structural tensor for the partic- 
ular case of a material with cubic symmetry and follow-up with a description 
of the polynomial basis of invariants constructed from this structural tensor 
and the right Cauchy-Green tensor C. 



3. 1 Structural tensor for cubic anisotropy 

To determine the structural tensors corresponding to crystals with cubic struc- 
ture we follow Zheng ([12]), where the structural tensors for all the different 
crystal classes are developed. Zheng's paper is based on the available prop- 
erties of Kronecker products of orthogonal transformations which allow the 
development of a simple method to determine the structural tensors with re- 
spect to any given symmetry group. As it has been highlighted in [12] there 
may exist many possible sets of the structural tensors for a given symmetry 
group, so one goal set by the author has been to find out the simplest irre- 
ducible representations; in particular, it is shown that each of the anisotropic 
symmetry groups can be characterized by a single structural tensor, and that 
is the result we will make use of. 

The crystal class corresponding to the Hexoctohedral cubic system is charac- 
terized by the following generators of its finite symmetry group 

Qi(|),Q2(|),Q3(|),-i 

where Q,i{6) refers to the rotation transformation about the axis e^, corre- 
sponding to a positive oriented orthonormal triad of vectors, through an angle 
6, and 1 is the second order identity tensor. 

Considering these group generators and applying properties of Kronecker 
products, the structural tensor for this particular cubic system is determined 
to be ([12]) 

M = ef +e^^^ +e'i^ (6) 

where the notation e^-^^ = Sj ® ® (8) e, has been introduced. The novelty of 
the approach suggested in this work is in the use of this fourth order structural 
tensor M in building up the energy function. All anisotropic structural tensors, 
investigated in the literature so far, have been second order. 
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3.2 Polynomial Basis 

An irreducible polynomial basis consists of a collection of members, where 
none of them can be expressed as a polynomial function of the others, i.e., 
they are independent scalars, and any other polynomial invariant of the same 
tensors can be written as a polynomial function of the basis members. The 
Hilbert theorem guarantees that a finite integrity basis exists for any finite 
basis of tensors ([5]). 

Taking [1] as a reference, we shall present an analogous procedure for the 
construction of specific constitutive equations based on functions whose argu- 
ments are the (joint) invariants of the right Cauchy-Green tensor C and the 
structural tensor M. Next we present the integrity basis invariants which will 
be the arguments of the constitutive functions to be proposed. 

The integrity basis consist of the traces of products of powers of the argument 
tensors. They can be divided in two main groups: the principal invariants, 
which involve invariants of the deformation tensor alone or the structural ten- 
sor alone, and the so-called mixed invariants, which consider joint invariants 
of both tensors. In the following we present separately the different kinds of 
invariants that can be formed from the right Cauchy-Green tensor C and the 
structural tensor M. 

• Invariants of the right Cauchy-Green tensor alone 

The principal invariants of the second order tensor C, denoted by = 
/fc(C), k = 1,2,3, are defined as the coefficients of the characteristic poly- 
nomial /(A) = det [Al — C] (see Appendix C for details). The explicit ex- 
pressions for the principal invariants of the second order tensor C are 

/i = tr(C) 
< /2 = tr(cofC) 
/3 = det (C) 

which can be expressed in terms of the basic invariants Ji, i = 1,2,3, 
defined as the traces of powers of C : 



Jl 


= tr(C) 


= 1 : C 




= ti{C^) 


= 1 : C2 






= 1 : C3 



• Mixed invariants of the right Cauchy-Green and the structural tensors 

In the case of several tensor variables, we use the term mixed invari- 
ant, even though the term simultaneous invariant can also be found in the 
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literature ([14]). 

We will follow Betten ([15]) to determine the scalar invariants of the 
tensors C and M. To construct a set of mixed invariants of the second-order 
tensor C and the fourth-order structural tensor M we consider a theorem 
presented in the reference [15] and its generalization for fourth-order tensors 
using the Hamilton-Cayley theorem, which means that powers of tensor 
and higher can be expressed in terms of 1, M,M^, M""-*-, where n 
represents the vector space dimension of C (for a symmetric second-order 
tensor n = 6). The additional mixed invariants are 

= C : M'^ : C 

ll = C:yi^ : C2 (7) 

Jg^ = : : C2, 

where k = 1,2,3,4,5. As shown in [15], the proof of the theorem relies on 
the assumption that M satisfies the symmetry conditions 

MijKL = MjiKL = MijLK = MklIJ. (8) 

Clearly these conditions are fulfilled by the structural tensor M = e^^^ + 

In addition to the invariants (7), [15] proved that the following expressions 
are also invariant 



Iii = l: M^- : C 
% = 1 : M^- : C^. 



(9) 



We will skip writing the superscript k in (7) and (9) when k = 1. 
Additionally we will also need the invariant 



1 : M : adj(C). 



Invariants of the fourth-order structural tensor alone 

The only remaining basic invariant of the single tensor M, formed as 
tr(M), is constant, and therefore it is not useful for the construction of 
strain energy functions. 



3.3 Poly convexity condition 



In this section we briefly describe sufficient (but not necessary) free energy 
function conditions which guarantee the existence of minimizers of some vari- 
ational principles for finite strains. As already mentioned, polyconvexity is the 
property of interest to us. 
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Local existence and uniqueness theorems in nonlinear elastostatics and elas- 
todynamics are based on strong ellipticity. The elhpticity condition states 
that the elastic free energy ip{F) leads to an elhptic system if and only if the 
Legendre-Hadamard condition 

VF,Ve,r/ e 3?' : /^i.^(F)(e®r/,e®r/) >0 

holds. 

The early global existence theory for elastostatics was based on convexity of 
the free energy function. However, that condition, as shown in [2], is unrea- 
sonable from a physical point of view. Using the notion of quasiconvexity due 
to Morrey ([3]), Ball ([2]) proved global existence theorems for elastostatics. 
In particular, it was proven that quasiconvexity implies the existence of min- 
imizers of some variational principles in finite elasticity. The quasiconvexity 
condition reads 

VF,Vcu G Co°°(S) V'(F)|S|= / ij(F)dV < [ ij(F + Vuj)dV 

Jb Jb 

Unfortunately this integral inequality is complicated to handle. A concept of 
greater practical importance is that of polyconvexity in the sense of Morrey 
([2]). Following Marsden ([13]), we say that the energy function ip is polyconvex 
if and only if there exits a function ip with arguments F, adj(F) = det(F)F~^ 
and det(F) such that 

^(F) = v.(F,adj(F),det(F)) (10) 
and if is convex function. 

As an illustrative example we present the case of ■?/'(F) = /(detF), for a convex 
function /. This function is not convex taken as a function of F (because the 
range of definition of F is not convex), however, it fulfills the polyconvex 
condition, since the condition of polyconvexity requires to take det(F) as the 
independent variable and, by hypothesis, the function / is convex in that 
variable. The polyconvexity condition has additive nature, i.e., if the functions 
ipi, i = 1,2,3 are all convex in their respective arguments then the function 
■ip(F) = ipi{F) + ip2{^d}(F)) + ?/'3(det(F)) is polyconvex. This property turns 
out to be very useful when proposing models because it permits to construct 
energy functions out of simpler ones. 

Due to the material indifference condition, the dependency on F of the energy 
function can be completely replaced by dependency on C, but the polycon- 
vexity condition does not translate to functions of C in a simple manner. 

Finally, we summarize the implication chain relating all the previous concepts 

convexity =^ polyconvexity ^ quasiconvexity =^ ellipticity 



11 



None of the opposite implications is true as counter-examples have been found 
for all of them ([13]). 



3.4 Isotropic free energy terms 

For completeness, here we present two statements about polyconvexity of some 
simple isotropic functions. The interested reader should consult [1] about fur- 
ther details on polyconvexity of various isotropic functions. 

Statement. The polynomial function 

F 7(tr(F^F))*^ = 7jf with A; > 1 and 7 > (11) 
is polyconvex. 

Proof. The function (tr(F^F))'= = ||F|p^ can be considered to be a function 
of F only and therefore, referring to the results in the previous section, it 
is enough to prove convexity relative to the argument F. As described in 
[1] the one possible approach to show convexity is to check the positivity of 
the second Gateaux derivative: 



<d(f^f)',h>=4 



F + eH)^ (F + eH) 
and from here the second derivative yields 

< (f^f)' ,iH,H)>=^<D (f^fY, H > 



e=0 



2k ||Ff^~^ < F,H >, 



e=0 



= 2k {WFf'^ < H, H > +{2k - 2) ||Ff '^"^ < F, H > 0. 

The desired result is established by noting that the constant 7 is positive 
and does not modify the signs of the derivatives. □ 

Statement. The functions 

F ^ 7(det(F^F))'= = 7/^, with k > 1, 7 > and 

F ^ -/51og(det(F^F)) = -/31og(/3), with /5 > 
are polyconvex. 

Proof. After noting that Is = (det(F))^, it is sufficient to show convexity rel- 
ative to det(F). The convexity is established by checking the non-negativity 
of the second derivatives of x'^'' and —2 log(x) which is a trivial exercise. □ 



12 



3.5 Anisotropic free energy terms 



Next we analyze the polyconvexity of some terms dependent on the structural 
tensor M. 

Statement. The polynomial function 

F ^ 7 (tr (f^FMF^f))'' = 7/4 with A; > 1 and 7 > 



is polyconvex. 

Proof. Mimicking the approach used for J^, we proceed by showing that /| 
is a convex function of F. Given that 



tr (F^FMF^F 



'F'^'F) : M : fF^F 



we can obtain the first and second Gateaux derivatives of /|: 



< D ((F^F) : M : (F^F) j , H > 

- A 

= 2k ((F^F) : M : (F^F))^"' (f^H : M : F^F + H^F : M : F^F 



;F + eH)'^ (F + eH) : M : (F + eH)^ (F + eH) 

fc-i 



e=0 



< ((F^F) : M : (F^F))'' , {H, H) >= < D ((F^F) : M : (F^F))\ H > 



fc-2 



F^H : M : F^F + H^F : M : F^F 



: Ak{k - 1) ((F^F) : M : (F^F) 
-2k ((F^F) : M : (F^F))''"^ x [2H^H : M : F^F + F^H : M : H^F- 
F^H : M : F^H + H^F : M : H^F + H^F : M : F^H] = 



= lQk{k - 1) ((F^F) : M : (F^F; 
+Ak ((F^F) : M : (F^F; 



fc-2 



F^H : M : F^F 



+ 



H^H : M : F^F + 2H^F : M : H^F 

■ (12) 

where the last equality used the symmetry properties of the structural tensor 
M, more specifically relations (B.l), (B.2) and (B.3). To complete the proof 
we separately analyze each term of the second derivative and show its non- 
negativity. 

The non-negativity of I4 follows from: 

3 

h = (F^F) : M : (F^F) = ^(F^F) : e^ (g) e^ : (F^F) = 

i=l 

3 

= 5^[F^F:e,®e,]2>0, 

i=l 



e=0 
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where property (B.4) has been used. In a similar manner, making use of 
both (B.4) and (B.5), we show that the other terms participating in the 
expression for the second derivative are also non-negative: 



H^H : M : F^F = ^(H^H : e,2)(e^ : F^F) = Y.{Heif{Feif > 0, 

i=l 1=1 
3 

H^F : M : H^F = ^(H^F : ® 6^)^ > 0. 



i=l 



The non-negativity of all terms in the expression for the second derivative 
together with the nonnegativity of 7 implies the desired result. □ 



Next, analogously to [1], we prove the polyconvexity of ^173 . This is the 
isochoric term corresponding to J4. 



r 

^3 



Statement. The function 



tr (f^FMF^f) I 



^ det(F)^/3 =^Jf with 7>0 



is polyconvex. 
Proof. Let </)(x, y) = and 



V',(F,C) = </>(||Fr/||,C)- 



^4/3 



where r/ is an arbitrary vector. We will establish that V',, is a convex function, 
when considered as a function of both arguments simultaneously. Condition 
(B.6) is satisfied for p = 4 and a = 4/3, therefore the function (j){x,y) is 
convex. The convexity of i/j^jiFX) follows from the following sequence of 
inequalities: 

m\FiV\\ , Ci) + (1 - m\\F2V\\ , C2) = AV',(Fi, Ci) + (1 - A)^A,(F2), 
where the triangular inequality and the convexity of y) have been used. 
The required result can be obtained directly from: 

h C : (ei^ + e2^ + e,^) : C ||Feif + \\Fe^\\' + HFeaf 



(det(F))4/3 (det(F))4/3 
= V^ei (F, det(F)) + V^ealF, det(F)) + ^,,{F, det(F)). 
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The positive coefficient 7 does not influence the conclusion. The same proof 
can be apphed to 7-^ as long as < a < 3. □ 

As pointed out in [1] the apparent symmetry between F and adj(F) in the 
definition (10) suggests that new polyconvex functions can be obtained by 
switching the deformation gradient tensor F with its adjoint tensor adj(F). 
The reader must be warned that replacing C with adj(C) is not equivalent to 
replacing F with adj(F) because 

adj(C) = det(C)C-i = det(F)2(F'^F)-i = det(F)2F-^F-'^ = adj(F)adj(F)^. 

The difference comes from the position of the transpose symbol: in C = F-^F 
the first matrix is transposed, while in adj(C) = adj(F)adj(F)-'" it is the second 
one, but the proofs already presented in this section are independent of the 
position of the transpose symbol and we will exploit this to verify the following 

Statement. Let = adj(C) : M : adj(C), then the functions 

Fh^7 (adj(F^F) : M : adj(F^F))'' = 77^3 with A; > 1 and 7 > 0, 

adj(F^F) : M : adj(F^F) K3 

F^l^ ,4/3 - = ^-^ with7>0 

-'3 -'s 

are polyconvex. 

Proof. If is a scalar invariant of C and M then the two functions will be 
polyconvex because the proofs of the previous two statements are indepen- 
dent of the position of the transpose symbol and will not be repeated here 
(for the second function, a = 8/3 is in the range of values which preserve the 

validity of the previous statement). The function is the proper isochoric 

^3 

variant of K3. 

We only need to prove that A'3 is a scalar invariant of C and M. We 
will prove on the way that Ki = adj(C) : M : C, K2 = adj(C) : M : 
and /m = 1 : M : adj(C) are also scalar invariants. If the characteristic 
polynomial (C.l) is multiplied as a double scalar product by C~^M : C, 
then the following equation for Ki is obtained 

Ki = h- hh + hlu. 

If Ki can be expressed as a combination of scalar invariants of M and C , 
then Ki is a scalar invariant itself. 

Proceeding in the same manner for A'2, one obtains after multiplication 
by C-^M : C2 

K2 = h — hh + hhi- 
For the multiplication is by C^^M : adj(C) and the result is 

K3 = K2- hK, + hlM. 
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Given that Ki and K2 are scalar invariants, for to be a scalar invariant 
one needs to show that Im is a scalar invariant. This can be accomplished 
by multiplication by C~^M : 1 with the result being 

Im = hi — hhi + h- 

In fact, for our specific tensor M = + 65 + 63 the invariant Im is equal 
to the invariant h- □ 



4 Model for polyconvex free energy function with cubic anisotropy 



Having presented proofs of polyconvexity for some strain energy functions, 
we proceed to propose a global model based on these functions for the case 
of materials with cubic anisotropy. The mathematical sanity of the model is 
guaranteed in advance by means of the polyconvexity of the proposed strain 
energy function. Consequently, this allows the application of the theorems 
concerning the existence of minimizing sequences discussed in Section 3.3. 

The proposed model is a linear combination of the energy functions proved 
above to be polyconvex. The general form of the free energy function reads 

V> = g (- log(/3) ) + + ih + Ih (13) 

where the last three terms have been selected to be as simple as possible, even 
though polyconvexity was shown for more general cases with exponents larger 
than 1. 

For completeness we briefly present the relationship between stresses and the 
free energy function based on the expression (4). Making use of the additively 
decoupled formulation of the polyconvex function (13) we expand the stress 
tensor in the following general manner 



From (14), following a procedure analogous to (3), we can deduce the formal 
expression for the tangent matrix: 

^ dS dCdS dS d"^^ 

C = — = = 2 — = 4 — - = 

dY dYdC dC dC^ 
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dij dC^ dljdh dCdc\ ^ ^ 

Expressions for each individual term corresponding to the first and second 
derivatives of the invariants can be found in Appendix A. 



4EE 

i=l 1=1 



5 Adittional Conditions 



In order to make a connection between the model and the physical data some 
conditions will be imposed on the model. These conditions will help determine 
uniquely the values of the arbitrary constants accompanying the individual 
functions appearing in the proposed model. 

The conditions discussed here refer to the comparison between the model re- 
sponse and the values of the parameters characterizing the physical behavior 
of a material in a natural state. The natural state chosen in this model corre- 
sponds to unstressed, undeformed configuration, i.e. C = 1. 

With the help of (14) and (15) the following two conditions can be formulated: 

• Stress free reference configuration 

This condition states that stresses must be zero when the deformation 
gradient becomes unity. Physically, this is equivalent to not having any 
remanent tensions when the material is totally unloaded. Mathematically, 
the stress free reference configuration means S(l) = 0, or upon substitution 
of the numerical values into (14) 

-a + /3 + 27 + 5 = 0. (16) 

• Tangent matrix at reference configuration 

The operation to be performed is the identification between the tangent 
matrix (15), particularized at the origin, C = 1, with the physical values 
corresponding to the classical elastic moduli matrix of a cubic material. 
There is an implicit assumption that the values of the elastic constants 
remain steady (time independent) for different values of the tensor C; albeit 
not being a completely realistic assumption, it can be considered a well- 
posed first approximation in the development of the model. Thus, we have 
the identification 



Co = 2dcS\^^. = Ad'c^p 



c=i 



lc=i ~ ^^c< 

where Cq represents the tangent matrix at the origin. Substitution of the 
numerical values for the derivatives gives the following three equations: 
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a + 27 = Cii (17) 
/? = Ci2 (18) 
a-p = 2Cu, (19) 

where Cn, C12 and C44 are the standard elasticity constants in Voigt notation. 
The solution of the system of equations (16-19) is: 

a = C12 + 2C44 

7= 2^^^^ ~ ^12 - 2C44) 
S = —Cii + C12 + 4C44. 

The nonnegativity of the elastic constants implies nonnegativity of a and (3. 
For 7 and 5 to be nonnegative the following condition must be satisfied: 

5 < < 1. (20) 

Conversely, if the inequalities (20) are satisfied then a, /3, 7 and 5 are non- 
negative. The ratio A = q'^^^q^^ is called the anisotropy ratio (see [16]) and its 
emergence in condition (20) suggests that the relevant measures of anisotropy 
are an integral part of the model itself. The data available in Hirth ([16]) 
shows that only five transitionary metals which are next to each other in the 
periodic table satisfy both inequalities: Cr, Mo, W, V, Nb. All five of them 
have body-centered cubic crystal structure. Some compound solids with cubic 
structure like AgBr and NaCl also satisfy the inequalities (20). 

The polyconvexity properties of /s and Iq are currently unknown to the au- 
thors, but the addition of terms proportional to or Jg in the strain energy 
function will not modify the condition (20) and consequently not enlarge the 
set of materials to which the model can be applied. 



6 Study in ID 



As an initial test for the model, two simple deformation gradient tensors have 
been applied as inputs: 



A 
1 
1 



and F 



1 7 
1 
1 
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stretch X Stretch X 



Fig. 1. Model response to simple stretch 




Shear y Shear y 



Fig. 2. Model response to simple shear 

Figure 1 displays the Cauchy stresses a as functions of the stretch A. Similarly 
to linear elasticity, the model predicts that the stresses an and (J22 = (^33 
are increasing functions of A. The apparent agreement between the proposed 
model and linear elasticity in the small strains regime is not surprising given 
that the model parameters are fitted to the hnear elasticity coefficients at zero 
strain. The physically desirable behavior of the stresses going to +00 (—00) 
as the stretch A goes to +00 (0) is also present. The behavior of the model in 
simple shear is shown in Figure 2. As in the previous figure, the stress curves 
predicted by the model are tangent to the stress curves for linear elasticity. 
The most notable difference in this case is that the model predicts nonzero 
stress (Jii while according to linear elasticity an = 0. 



7 Variational formulation and finite element discretization 

Consider a body B with boundary dB = AiU A2. The boundary Ai consists 
of all surface points where displacement is applied and the boundary A2 of all 
surface points where tractions are applied (.4i nA2 = (p)- The boundary value 
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problem can be formulated as (following [17]): 



Div(P) = in;B, (21) 

P = P(Vu,x) in;B, (22) 

u = u on (23) 

Pn = t on ^2, (24) 



where (21) expresses the equilibrium condition of the solid in the absence 
of body forces, (22) is the constitutive law for the solid, (23) expresses the 
boundary condition for the section of the boundary Ai on which displace- 
ment is imposed and (24) refers to the section of the boundary A2 on which 
the tractions are imposed. For hyperelastic materials possessing strain energy 
function ijj the first Piola-Kirchhoff stress tensor is given by P = II where 
F = Vu. 

It can be shown ([17]) that the solution u of the problem posed by (21-24) 
in the case of polyconvex strain energy function ip is the minimizer of the 
functional 

J(u) = / ^(x, u, \/u)dV - [ u ■ tdS, (25) 

Jb JdA2 

which can be used to formulate the finite element discretization. The spatial 
discretization is accomplished by representing the body B as a. union of disjoint 
elements, B = U^'f" ^e- Even though many different elements can be used 
for the discretization, for simplicity, we will assume that the discretization has 
been achieved by the use of second order tetrahedral elements. The unknowns 
to be solved for are the nodal displacements For each element, the dis- 
placement is expressed as a sum of the nodal functions, u^^^^ = J2]!Li ^jf^'^'^a'^- 
The global displacement approximation becomes = Y^^=i"^ Y^^^iN^^^u'^^K 
Substitution of into (25) leads to a discrete functional Jh{uh). The mini- 
mization procedure for the discrete functional gives rise to a system of equa- 
tions which can be solved for the nodal displacement. Further details about 
the finite element procedure can be found in [18]. 



8 Numerical examples 



In this section we will consider two basic examples which illustrate the agree- 
ment of the model with linear theory at small strains and the departure from 
it at large deformations. 
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Fig. 3. Stress concentration and convergence plots for plate with an initially circular 
hole 

8.1 2D Example - Plate with Hole 



To verify the consistency of our model with linear theory and to check its 
convergence we considered the problem of uniaxial stress applied to a plate 
with initially circular hole. In addition to the well known analytical solu- 
tion for isotropic linear elastic material, this problem has analytic solution 
for orthotropic linear elastic materials [19]. The stress concentration factor, 
calculated from this solution specialized to cubic anisotropy with one symme- 
try axis perpendicular to the plane of the plate and another symmetry axis 
aligned with the loading direction, is shown on the left side of Figure 3. The 
numerical solution obtained from our model is in good agreement with this 
analytical solution when the applied stress is small compared to the elastic 
constants of the material. As the stress is increased to become comparable to 
the elastic constants the nonlinearities become important and approximately 
the same reduction in the minimum and the maximum stress concentration is 
observed. 

The plot on the right of Figure 3 shows the dependence of the error on the 
mesh size. Within the linear regime the convergence is quadratic as expected 
for Newton-Raphson solvers. The convergence is somewhat reduced for the 
nonlinear regime, but it is still within acceptable levels. 



8.2 3D Example - Circular Bar 



The problem which we will consider here is the extension of a single crystal 
cylindrical bar. The bottom end of the bar is clamped and the side surface of 
the bar is traction free. The top end is displaced in the axial direction by a 
specified amount, but it is left free to move in the plane perpendicular to the 
original bar axis. Three different orientations of the bar axis relative to the 
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crystal will be considered. In each case the bar axis will coincide with one of 
the following crystallographic directions: [100], [110] and [122]. Top and side 
views of the deformed bar for the three cases are shown in Figure 4. For better 
visualization the applied displacement is equal to 100% of the bar length, but 
the behavior is similar at smaller stretches. 

In the first case the cross-section of the bar remains approximately circular 
and the extension process is approximately symmetric about the bar axis. 
Similar behavior is observed if the bar axis is aligned with the [111] direction. 

The anisotropic response is clearly observed in the second case. The contrac- 
tion of cross-section of the bar in the two directions is markedly different due 
to the anisotropy brought by the term containing I4. As it could be expected 
the contraction in the [001] direction is much less than in the [110] direction. 
A possible interpretation of this effect for metallic lattices is that part of the 
contraction in the latter direction is accomplished by atomic bond rotation 
rather than pure extension/contraction of the bonds. 

The third case illustrates the development of transverse displacements when 
the crystal lattice lacks enough symmetries relative to the loading axis. The 
tilting effect is completely due to the anisotropy. If the movement of the top 
end is restricted significant transverse stresses will develop. 



9 Conclusions 



A new model for materials with cubic anisotropy has been proposed in this 
paper. The model is based on additively decoupled strain energy function 
which satisfy the polyconvexity condition and therefore guarantees existence 
of minimizing sequences for the appropriate variational functionals. The poly- 
convexity of new strain energy terms capturing the anisotropy of cubically 
symmetric systems has been shown. A simple strain energy function capable 
of capturing the fundamental effect of the anisotropy ratio has been suggested 
and tested in numerical simulations which reveal that the model possesses 
many of the relevant physically desirable properties. 

The main difference between this model and the orthotropic models in the 
literature (for example, [1]) is the use of a single fourth order structural tensor. 
In spite of the complications coming from the higher order of the tensor, 
the model avoids a major complication which the orthotropic models face: 
enforcing the equality of the properties in the three mutually perpendicular 
symmetry directions. 
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A Invariants derivatives 



The computation of the stress tensor S and the tangent tensor C requires the 
first and second derivatives of the invariants. Expressions for these derivatives 
deduced in the reference configuration are given below. 



• First derivatives 



dh 
dCij 

dh 
dCu 

dh 
dCjj 

dh _1 

dh _1 
dCjj 2 



OlJ 

h5ij - Cij 

{MjJpQ + MJIPQ + MpQJJ + MpQJj) CpQ 

(SimSjn + SinSjm) MmnpqCpq+ 



CmnMmnpq {SipCjQ + Ciq5jp + Cpj5qj + SiqCjp)] 
1 

dCij 2 



dh 1 I 

{SimCjn + SjmCin + Cim6jn + Cjm^in) MmnpqC pq+ 



CmnMmnPQ {SjpCjQ + CjqSjp + CpiSqj + SiqCjp) 

dhi 1 

— = -^mnMmnpq {Spi^Qj + Spj^Qi) 

OUij I 

dhi 1 

-Kp^ — = -^MnMmnPQ {SpiCqj + SpjCgi + CpjSqj + CpjSqi) 
OUij z 

= ^hhmMMNPQ {^CpqCj} - Cp]Cq] - CpjCg]) 
{Mkljp + Mlkjp) Cpi + {Mklpj + Mlkpj) Cpi + 

CpQ {MpQiiSjK + MpqjkSjl) + 



• Second derivatives 
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dCijdCKL 
dCudCKL 
dCijdCKL 



dCijdCKL 4 



dCjjdCKL 



dCjjdCKL 



M 



dCjjdCKL 

hi 



dCijdCKL 



M 



JK 



1 

[MijKL + MjiKL + MijLK + MjiLK + 

Mklij + Mklji + Mlkij + Mlkji] 

■- \ [{MijKP + MjiKp) CpL + {MijPK + Mjipk) CpL + 
{MijpL + Mjipl) CpK + {MijLP + MjiLp) CpK + 

{Mklip + Mlkip) Cpj + (MxLP/ + Mlkpi) Cpj + 
(Mt^ljp + Mlkjp) Cpi + (M/^LPj + Mlkpj) Cpi + 

CpQ {Mpqil5.jk + MpqjkSjl) + CpQ {MpQiidjK + MpqjkSjl) + 
CpQ {Mpqil5.jk + MpqjkSjl) + CpQ (MpQKlSjL + MpqliSjk)] 

■- ^ [(M/ipQ^j/^ + MjkpqSjl + MjlpqSjk + MjkpqSil + 

+MkjpqSilMljpq5ik + Mkipq5jl + MlipqSjk) CpQ + 
Clp {MrpiqCjq + MxpjqCqi + MkpqjCqi + MkpqiCq 

ClP {MpkIqCjq + MpkjqCqi + MpKQjCqi + MpkqiCq 
CrP {MlpjqCjq + MlpjqCqi + MlpqjCqi + MlpqjCqj 
CkP {MpliqCjq + MpljqCqi + MpLQjCqi + MplqjCqj) + 
CpQ {MpqilSjk + MpqjkSjl + MpQjiSjK + MpQKjSiL + 
MpQLjSiK + MpQKlSjL + MpQLlSjK + MpqjkSil)] 



QJ) 

'qj) 

Qj) - 



-SmnMmNPQ [Spi (SqkSjL + SqlSjk) + Sqi (SpK^jL + SplSjk] 
5pj (SikSql + SilSqk) + ^QJ (SikSpl + ^/L^Pi^)] 

1 



Iz^mnMmnpq 



^■-^/J '-'PQ ^ ^PI^QJ " ^PJ^QI ) ^KL' 



PQ l^'-^/x'-' JL + ^IL ^JK ) "~ '-'/J yyPK^QL + ^PL^QK 



1 1 
1 1 

^'-^PJ XyQK^IL "T ^QL^IKJ "T 2 \^'-"PA''-^,/L "r ^PL^JK 
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B Proofs of basic properties 



Statement. 

H^F : M : F^H = H^F : M : H^F (B.l) 
Proof. Making use of the indicial notation we liave 

H^F : M : F^H = HmiFmjMjjklFnkHnl- 

Considering the symmetries of M, expression (8), and redefining properly 
dummy indices, the above expression can be rewritten as 

H^F : M : F^H = HmiFmjMulkFnrHnl = HmiFmjMulkHnlFnk 
= H^F : M : H^F.D 

Statement. 

F^H : M : H^F = H^F : M : H^F (B.2) 
Proof. Analogously to the previous statement, we have 

F^H : M : H-^F = FmiHmjMjjklHnkFnl = FmiHmjMjjklFnrHnl 
= HmjFmiMjiklHnrFnl = H^F : M : H^F.D 

Statement. 

H^H : M : F^F = F^F : M : H^H (B.3) 
Proof. Similarly to the two statements above, we have 

H^H : M : F'^F = HmiHmjMjjklFnrFnl = FnrFnlMjjkiHmiHmj 
= FnkFnlMklijHmrHml = F^F : M : H^H.D 

Statement. If e is a vector and A, B and C are second order tensors then 

A : B®B : C = (A : B)(B : C) and (B.4) 
e® e : A^A = (Ae)2. (B.5) 

Proof. Using index notation, we have the following equalities: 

A : B ® B : C = AjjBjjBklCrl = (A : B)(B : C) 

and 

e ® e : A^A = ciCjAkiAkj = AkiCiAkjCj = {Aef.n 
Statement. If a > and p >2 satisfy the condition 

> B.6 

a p — 1 
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then the function </)(x, ?/) = ^ is convex when considered as a function of 
both arguments simultaneously. 
Proof. The complete proof is given in [1] . We will only note here that the con- 
dition (B.6) is equivalent to the positive definiteness of the second derivative 
(the Hessian) of 0. □ 



C Characteristic polynomial 

The characteristic polynomial, also known as Cayley-Hamilton polynomial, 
for a matrix C G M^^^, is the equation resulting from the eigenvalue problem 
corresponding to that matrix C, i.e., det(Al — C) = 0, giving 

= - tr(C)A2 + tr(adj(C))A - det(C), 

or equivalently 

= - tr(C)C^ + tr(adj(C))C - det(C)l. (C.l) 
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